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ABSTRACT 

In this paper we discuss the applications of high order compact finite difference methods 
for shock calculations. The main idea is the definition of a local mean which serves as a 
reference for introducing a local nonlinear limiting to control spurious numerical oscillations 
while keeping the formal accuracy of the scheme. For scalar conservation laws, the resulting 
schemes can be proven total variation stable in one space dimension and maximum norm 
stable in multiple space dimensions. Numerical examples are shown to verify accuracy and 
stability of such schemes for problems containing shocks. The idea in this paper can also be 
applied to other implicit schemes such as the continuous Galerkin finite element methods. 
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1 Introduction 


Compact schemes are methods where the derivatives are approximated not by polynomial 
operators but by rational function operators on the discrete solutions. In this paper we are 
interested in solving a hyperbolic conservation law 


u t + f(u) x + g(u) y = 0 

u(x,y,0) = u°(x,y) (1.1) 

using compact schemes. In the semi-discrete form, a compact scheme for solving (1.1) can 
be written as 

u))„ - (1-2) 

where A and B are both local , one dimensional operators. The subscript x or y indicates 
that the operator is applied in the x or y direction. 

For example, a fourth order central compact scheme is given by (1.2) with 

(v4v), = ^( u i-i + + u «+i) 

(Bv)i = i(v,- + i - u,-i), (1-3) 

a sixth order central compact scheme is given by 

(Au)i = -(u«-i T 3 U « + u «'+i) 

( Bv)i = -^(vi +2 + 28v i+ i -28vi-i - Vi- 2 ), (1-4) 

and two third order upwind compact schemes are given by 

(Au)i = ^( _U '-1 + 5u i - U *’+l) 

(Bv)i = ^(3u< - + u,_ 2 ) (1.5) 

and 

(Au), = h-Vi-i + ~ u *+i) 

( Bv)i = -(— Uj + 2 + 4u,+i — 3u t ) (1*6) 

depending upon the wind direction. Notice that (1.5) and (1.6) have the same implicit part 
A which is symmetric. This fact will be used later in Section 2 to define our local means. 
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The cost of compact schemes, regardless of the number of space dimensions, involves only 
inversion of the narrowly banded (usually tridiagonal) matrix A and hence is comparable 
to explicit methods. This is notably different from other implicit methods such as the 
continuous Galerkin finite element methods in multiple space dimensions, even if they are 
similar in one space dimension. 

The advantages of compact schemes include the relatively high order of accuracy using a 
compact stencil (for example, the fourth order scheme (1.3) when discretized in time using 
Euler forward, uses only a three point stencil in each time level), a better (linear) stability, 
and usually fewer boundary points to handle. In recent years compact schemes have attracted 
much attention in various fields such as the direct numerical simulations of turbulence. We 
refer the readers to [2], [3], [4], [12], and [18] for more details. 

The objective of this paper is to apply compact schemes for shock calculations. As 
with any other linear schemes (schemes which are linear when applied to linear equations), 
compact schemes usually demonstrate nonlinear instability when applied to discontinuous 
data. We follow the TVD (total variation diminishing) ideas in [9], [13] and try to define 
a suitable nonlinear local limiting to avoid spurious oscillations while keeping the formal 
accuracy of the scheme. Notice that the compact scheme, like an implicit scheme, is global. 
That is, the approximation to f{u) x at x = X{ involves Uk along the whole line due to the 
tridiagonal inversion A~ x . Our main idea is to define a local mean, and use it as a reference for 
introducing a local limiting. In Section 2 we introduce the limiting for one space dimension 
and prove total variation stability. In Section 3 we introduce the limiting for multiple space 
dimensions and prove maximum norm stability. In Section 4 we present numerical examples, 
and concluding remarks are included in Section 5. 

The ideas in this paper were first used by us for continuous Galerkin finite element method 
in [7]. That is an on-going project. In this paper we restrict our attention to scalar problems 
in order to obtain provable stability results. The application of the method to systems of 
hyperbolic conservation laws and to other types of compact schemes (e.g. [1]) is currently 
under investigation. 

In this paper, we use the total variation diminishing (TVD) Runge-Kutta type time dis- 
cretization, introduced in [14], [17], to discretize the ODE in the method-of-lines formulation 
(1.2). In the second order case, the time discretization is 

«<»> = u n + AtL(u n ) 

«" +1 = \u n + l -uU + ^AtL(uW), (1.7) 

and in the third order case it is 

u (i) = u n + AtL(u n ) 

= ju” -|- + \AtL(u^) (1.8) 

4 4 4 

“” +I = ju" + f« (2 > + |a(Z( u < !) ). 

These special Runge-Kutta type time discretizations are labelled TVD because it can be 
proven that under suitable restrictions on the time step At (the CFL condition), the full 
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discretization (1.7) or (1.8) is TVD, or stable under another norm (for example, the L 0 0 
norm) if the first order Euler forward time discretization for (1.2) 


u n+1 =u n + A tL(u n ) 

(1.9) 

is TVD or stable under the other norm. For details, see [14] and [17]. 

We thus only need to consider the Euler forward scheme (1.9) for stability analysis in the 
subsequent sections. 

2 One Space Dimension 


In one space dimension, equation (1.1) becomes 


Ut “I - y(n)a : — 0 

u(x,0) = u°(x), 

(2.1) 

the scheme (1.2) is 



(2.2) 

and the Euler forward time discretization (1.9) becomes 


u" +1 = < + A tL{u n )i. 

(2.3) 

Scheme (2.3) can be easily written into a conservation form 


suitable for shock calculations. However, the numerical flux h n , i 
due to the tridiagonal inversion A ~ l . If we define 

(2.4) 

is not a local function of 

Ui = ( Au)i , 

(2.5) 

then scheme (2.3) can be left-multiplied by A to become 


< +l = «? - 

(2.6) 

When written into a conservation form, 


r-.n+l -n ( rn fn \ 

(2.7) 

this involves a numerical flux f? \ which is a local function of u£. 

For example, 

f i+ i = ^(/(«i+i) + /(«,-)) 

(2.8) 


for the fourth order central scheme (1.3), 
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(2.9) 


/,+i = ^ (/(«<+ 2 ) + 29/(u, + 1 ) + 29/(ti f ) + /(«i-0) 
for the sixth order central scheme (1.4), and 


3/M - /(<*-i)) 

(2.10) 

and 


/«+§• = 2 ( _ /( U *'+ 2 ) + 3/(u, + l)) 

(2.11) 

for the two third order upwind schemes (1.5) and (1.6), respectively. Notice that scheme 
(2.7) resembles a cell-averaged (finite volume) scheme [11], The u; in (2.5), like a cell average, 
is a local mean of u, defined by Au in (1.3) through (1.6). Since the computation of the flux 
f i+ 1 in (2.7) involves the values of u, a “reconstruction” from u to u 

Ui = {A * 1 *), 

(2.12) 

is needed. This reconstruction is global. 

It is now rather straightforward to define the limiting. We first write 


m = / + (t.) + /-( u> 

(2.13) 

with the requirement that 


a/ . +( “) >0 , a C (w) < 0 . 

ou ou 

(2.14) 

The purpose of this flux splitting is for easier upwinding at later stages. The simplest such 
splitting is due to Lax-Friedrichs 

/*(«) = \ (/(«) ±H> « = max \f'(u)\ 

(2.15) 


where the maximum is taken over the range of u°(x). We then write the flux f i+ i in (2.7) 
also as 


A A A 

/,•+! = + /;; 


where /* , are obtained by putting superscripts ± in (2.8) through (2.11). 

H" 2 

Next we define 


(2.16) 


= /itj - df- +i = /“(«,+,) - (2.17) 

Here df t are the differences between the numerical fluxes /* , and the first order, upwind 

fluxes / + (u,) and /“(fq+i). These differences are subject to limiting for nonlinear stability. 
We define the limiting by 
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= m(i/+L,A + / + (S j ),A + / + (fli-.)) 
dftp = m(^;,,A + /-(fl,),A + r(“i + i)) (2-18) 

where A+i>,- = v,-+i — u,- is the usual forward difference operator, and the (now standard) 
minmod function m is defined by 


m(aj, •••,<**) = 


smini<,<*. |a,|, if sign(ai) - ■ • • = sign(a k ) = s 
0, otherwise. 


Notice that the limiting defined in (2.18) is upwind biased. 
The limited numerical fluxes are then defined by 


= / + (8 ( ) + df* /-'f' = /-(«.+.) - df -'J 


P+(m). 


r~( m ) _ 


?-(m) 


and 


(2.19) 


( 2 . 20 ) 


f( m ) _ f+( m ) i r - ( m ) 

A'+J “ •'i+i + J i+\ ■ 

If we define the total variation of the mean u by 


(2.21) 


TV(S) = £>*, - 5.-I (2-22) 

% 

we have the following proposition. 

Proposition 2.1 

Scheme (2.7) with the flux (2.21) is TVDM (total variation diminishing in the means) 


TV(u n+l ) < TV(u n ) 

(2.23) 

under the CFL condition 


max (f + '(u) — f~'{u)) < jL 

mini u"<u<maxi u" ' A(C 2 

(2.24) 

Proof: We follow Harten [9] and write the flux difference as 


- ftl = -Cj+iA+a,- + A-i A + fi,_, 

(2.25) 

where 



c >*\ 


D 


*“2 


A + /-(« i )-^; { r ) +c l r’ 

2 * 2 

A + U; 

a +/+(«._ i) + d/jr 1 - df*} r° 

2 * 2 

A + u,_i 


(2.26) 
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The limiting in (2.18) and the properties of f ± (u) in (2.14) clearly imply 


C, + t > 0, D,_ r > 0 (2.27) 

and 

At (r IP w At /-2A + /-(ii,) + 2A + /+(ii j ) 
aJ v c ' + i + A+iJ s ^ 5^; 

The last inequality is due to the CFL condition (2.24). TVDM (2.23) is now immediate 
according to Harten [9]. 


J < 1 (2.28) 


□ 


In order to obtain total variation stability for u, we need the following simple lemma. 

Lemma 2.2 

If there are two numbers 0 < 8 < 1 and a > 0, which are independent of N, such that 
the N x N matrix A = (a,j) satisfies: 


1 N 

i<j<N IH “ and ■ l a bl < 6 \ a n \ > i = 1, • • ■ , iv 

n i = 1 

* 4- j 


(2.29) 


(strongly diagonally dominance for the transpose of A), then the L\ norm of A -1 is bounded 
independently of N, 


IU, < 


a 


1 -S' 


Proof : Let A = diag(a n , • • • , o^n), B — A — A and C = —BA -1 . We have 


(2.30) 


N 


\C\\l x = max ^ c,j = max V) - — — < 8 . 
1 Jl i<i<N \ ajj \ _ 

1=1 

* j 


Hence it follows that 


P _1 lk, 


< 

< 


||[(/-C')A]- 1 || Li =||A- 1 (/-C)- 1 || Li 

HA-MIl, Il(/-C)- 1 11L 1 <11A- I ik l _ [ | q[ ^ 

a 

1 - 8' 


□ 


6 



For most compact methods, the matrix A satisfies the condition (2.29) for Lemma 2.2. 
For example, in the schemes defined by (1.3), (1.4), (1.5) and (1.6), A satisfies the condition 
(2.29) with <5 = |, a = 6; 6 = §, a = 5; 5 = |, a = 3 and 8 = f, a = 3, respectively. For 
such compact schemes, we can now prove the total variation stability for u. 

Proposition 2.3 

If a compact scheme (2.7) satisfies the conditions in Proposition 2.1 and Lemma 2.2, then 
it is TVB (total variation bounded). That is, 

TV(u") = Z |< +1 - <| < C (2.31) 

for all n > 0 and A t > 0. Here C is a constant independent of n and A t. 

Proof : 

By (2.12), we have 

TV(u") = = 

x i 

< p-'lk E NVi - «?| < -rri TV (v\ 

% 

□ 

This Proposition guarantees convergence of at least a subsequence of the numerical so- 
lution. 

We now discuss whether the limiting defined in (2.18) maintains the formal accuracy 
of the compact schemes in smooth regions of the solution. For this we need the following 
assumption. 

Assumption 2.4 

Hi = ( Au)i = m + 0(Ax 2 ) (2.32) 

for all u 6 C 2 . 

□ 

This Assumption is satisfied by any compact scheme with a symmetric A , for example 
all those listed in (1.3) through (1.6). 

Under Assumption 2.4, it is easy to verify by simple Taylor expansions that 


A + / ± (ufc) = f ± (ui) x Ax + 0(Ax 2 ) k = i — l,i,i + 1 

dftt = ^/ ± (u,) I Aa: + 0(Ax 2 ). (2.33) 

Hence in smooth regions away from critical points (critical points are defined here as 
points for which f + (u) x = 0 or / - (u) x = 0), the second and third arguments of the minmod 
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functions in (2.18) are asymptotically of the same sign as the first argument and half in 
magnitude. Hence the first argument will be picked by the minmod function (2.19) for 
sufficiently small Ax, thus yielding 


(2.34) 

This guarantees the original high order accuracy of the scheme in such smooth, monotone 
regions. At critical points, the accuracy will degenerate to first order as a generic restriction 
of all TVD schemes (see, for example, [13]). To overcome this difficulty, we use a modification 
of the minmod function 




aj, if [cii | < M Ax 2 

m(ai, • • • , a,k), otherwise 


(2.35) 


where M is a constant independent of Ax. This modification is discussed in detail in [5] and 

[15]- 

With this modification we can obtain schemes which are formally of uniform high order 
accuracy and equal the original unlimited scheme in smooth regions including local extrema. 
Moreover, we can prove the following proposition. 


Proposition 2.5 

The conclusions of Proposition 2.1 and 2.3 are still valid for any n and At such that 
0 < nAt < T, with TVDM in (2.23) replaced by TVBM (total variation bounded in the 
means) 


TV(u n ) < C (2.36) 

where C is independent of At, if the minmod function m in (2.18) is replaced by the modified 
minmod function fh defined in (2.35). 

Proof: 

The proof is similar to that contained in [15] and [5] and is thus omitted. 

□ 

The choice of the constant M in (2.35) is related to the second derivative of the solution 
near smooth extrema. For details, see [5] and [15]. The numerical result is usually not 
sensitive to the variation of M in a large range. 

In this paper we only consider pure initial value problems. u° in (1.1) is assumed to 
be either periodic or compactly supported. For initial boundary value problems, u in (2.5) 
is defined differently at the boundary, as is the scheme (2.6). The limiting (2.18) can be 
modified at the boundary so that the scheme remains TVDM (or TVBM) and TVB for 
initial boundary value problems. We refer the readers to [5] and [16] for more details. 


3 Multiple Space Dimensions 

For notational simplicity we only consider the two dimensional case (1 . 1 )- (1 .2) . Three space 
dimensions do not pose additional conceptional difficulties. As before, we only need to 
consider the Euler forward time discretization 
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(3.1) 


utf' = «5 + AtL(u n )ij. 

We again define 


Uij = ( A y A x u)ij (3.2) 

so that scheme (3.1) can be left-multiplied by A y A x to become 

= < ~ ( 3 - 3 > 

Here and in what follows we will use the commutativity of A Xy A yy B x and B y so that a 
product can be written in any order. Scheme (3.3) can be written into a conservation form 


u 


"+ 1 = u n - —(f n 1 . - f n 1 .) - — (o" 1 - 0 n . i) 

•3 U,J AX^ J,+ *’ J AlA^+2 y, ’3~2 J 


Ay 


(3.4) 


which involves numerical fluxes f n i . and g n . , i as local functions of u^. For example, 

M 2 W bjT 2 


fi+y = (/(ui+ij) + /(«<,-)) 

5<J+| = (p(«i,j+i) + 5(“»j)) ( 3 - 5 ) 

for the fourth order central scheme (1.3), with analogous definitions for the other schemes. 
Again, scheme (3.4) resembles a cell-averaged (finite volume) scheme [10]. The u,y defined 
by (3.2) is a local mean of u, and a “reconstruction” from u to u 

Uij = {K' A y X u)ij (3-6) 

is needed to compute the fluxes /,+ ij and in (3.4). 

We remark that the additional costs of implementing scheme (3.4), comparing with the 
original scheme (3.1), are the two local operators A x and A y . The major part of the cost still 
consists of the two tridiagonal inversions. 

The limiting to obtain nonlinear stability can now be defined in a dimension by dimension 
fashion; we can use the one-dimensional flux splitting (2.13), for f(u), to write the flux /,+ij 
as 



+ f i+ i 

•i o >. 


(3.7) 


where , . are again obtained by putting superscripts ± in, for example (3.5). The remain- 
ing definition of the limiting parallels that in Section 2, with a dummy index j added for 
the reference y value. We still start with the differences between the high order numerical 
fluxes and the first order upwind fluxes 


and limit them by 


dyj = /'(«.'+u) - /r+ij 


(3.8) 
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(3.9) 


d %l) = "> (tfcy.W (s.' + w)) 

where = v i+i,j — is the forward difference operator in the x direction and the 

minmod function m is defined by (2.19). We then obtain the limited numerical fluxes by 


J i+y ~ } wu) + d / f+ ip 


i .+1 !! ~ f ("•+■■)) 




(3.10) 


and 


m ) _ f+( m ) I r-(n» 

+ j.i 


The flux in the y-direction is defined analogously. 

In light of [8] this scheme cannot be TVD in two space dimensions, 
obtain maximum norm stability through the following proposition. 


(3.11) 


However we can 


Proposition 3.1 

Scheme (3.4) with the flux (3.11) satisfies a maximum principle in the means 


under the CFL condition 


max u” +1 < max u" 
*<} J i,j 


(3.12) 


[max (/ + '(u)) + max (-/ '(«))] ~ + [max (g*‘(u)) + max (-<r'(tij)] < I (3.13) 

where the maximum is taken in min tJ - <u< max,^ u^. 

Proof : Similar to the development in Proposition 2.1, we can write the flux differences as 


& i + A,--i 


with 


A+Jd > 0, D t .y > 0, C iji+i > 0, 
due to the flux splitting (3.7), the limiting (3.9), and 


Au-J s 0 


Ax { C i+l: + D i-y) + ( C i,i+L + - 1 


the CFL condition (3.13). 
We then have 


(3.14) 


(3.15) 


(3.16) 
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«" +1 

u t; 


= oj + #; (A^AX; - (c jJ+i A» + aj - Aj-iAW,.,) 


1 - 


Ax 

At 


Ay 


u 


*3 


At 

Ax V' ,+ ?’ J ' "•-»•»/ Ay 

+ ^tj+prj+i + Aj-pij-i 




which implies the maximum principle (3.12) because u?* 1 is written as a convex combination 
of u? , u" ± jj and u" J±1 with positive coefficients which add up to one. 


□ 

In order to obtain maximum norm stability for u , we need a lemma similar to Lemma 2.2. 

Lemma 3.2 

If there are two numbers 0 < 6 < 1 and a > 0, which are independent of N, such that 
the N x N matrix A = (a,ij) satisfies: 

1 N 

max t — r < a, and V) la;,- 1 < 8 lad , i = (3.17) 

i<i<N |a,,| - “ 

j ¥■ * 

(strongly diagonally dominance for A), then the Loo norm of A -1 is bounded independent 

of N 


ll^-Mkoo < — 6 - 

Proof : The proof is similar to that for Lemma 2.2 and is thus omitted. 


(3.18) 


□ 

For the compact methods we consider, the matrix A is symmetric. Hence the requirements 
(2.29) and (3.17) are the same. 

We can now use Lemma 3.2 to obtain the maximum norm stability for u. 

Proposition 3.3 

If a compact scheme (3.4) satisfies the conditions in Proposition 3.1 and Lemma 3.2 for 
both A x and A v , then it is stable in the maximum norm. That is, 

maxli&l < C (3.19) 

»,j I J I 

for all n > 0 and At > 0. Here C is a constant independent of n and At. 

Proof : 

By (3.2), we have 
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max U” 


5 (tT7)V|°«I- 


U 


□ 

This Proposition does not guarantee convergence, but it at least guarantees that the 
numerical solution will not blow up due to instability. 

Under the Assumption 2.4 for both A x and A y , we can again easily verify that the 
limiting (3.9) maintains formally the original high order accuracy of the scheme in smooth, 
monotone regions. The degeneracy of accuracy at critical points can once again be overcome 
by adopting the modified minmod function (2.35) in the limiting (3.9). 


4 Numerical Examples 

To test the behavior of the schemes discussed in Sections 2 and 3, we use the one and two 
dimensional Burgers equation with the smooth initial conditions 


0 

0.3 + 0.7 sin(a;) (4.1) 



u(x,y,0) = 0.3 + 0.7sin(a: + y). (4.2) 

Both are assumed to have 27r-periodic boundary conditions. The solutions will stay smooth 
initially, and then develop shocks which move with time. The exact solution to (4.1) can be 
obtained by following the characteristics and solving the resulting nonlinear equation using 
Newton iteration. The exact solution to (4.2) is that of (4.1) with x replaced by x -+- y and 
t replaced by 2 1. These are standard test problems for scalar nonlinear conservation laws 
containing shocks. For comparison with finite difference ENO schemes and finite element 
discontinuous Galerkin methods, see [17], [5] and [6], 

The schemes we test are based on the fourth order central scheme (1.3) coupled with a 
fourth order Runge-Kutta time discretization (henceforth referred to as the central scheme), 
as well as the third order upwind schemes (1.5)-(1.6) coupled with the third order TVD 
Runge-Kutta time discretization (1.8) (henceforth referred to as the upwind scheme). For 
the flux splitting (2.13) we use the Lax-Friedrichs splitting (2.15). The time step At" is 
taken to satisfy a CFL condition 



u(x,0) = 


and 
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in one dimension and 


, , A t n 

max u? — — < 0.5 
i Ax ~ 


(4.3) 


max 



A t n A t n \ 

Ax ' Ay ) 


<0.5 


(4.4) 


in two dimensions. When the modified minmod limiter (2.35) is used, the constant M is 
taken as 1. 

We first test the effect of limiters when the solution is smooth but not monotone. In 
Figure 1 we plot the L\ error versus number of grid points, in a log-log scale, at t = 0.6 for 
the one dimensional case and at t = 0.3 for the two dimensional case. In such scales, the error 
should be a straight line with slope —k for a A:-th order method. We can see that the original 
compact schemes and the schemes with modified minmod limiter (2.35) (henceforth referred 
to as the TVB limiter) give the expected third and fourth order accuracy respectively, while 
the schemes with the minmod limiter (2.19) (henceforth referred to as the TVD limiter) give 
only second order accuracy due to the degeneracy at the critical points. We can also see 
that both the central and the upwind schemes work well for this smooth problem. 

We then test the effect of limiters when the solution becomes discontinuous. In Figure 2 
we show the results of the original compact schemes at t = 2 for the one dimensional case. 
We can see over- and under-shoots as well as oscillations, and in this case the result of the 
central scheme is much worse than that of the upwind one. In Figures 3 and 4 we show 
the results with the TVD and the TVB limiters. Apparently the limiters have stabilized 
the solution, as predicted by the theory. However the result with the central scheme is not 
quite satisfactory. In Figures 5 and 6, we show the pointwise errors, in a logarithm scale, for 
the numbers of grid points N = 10,20,40,80 and 160. We can see that the central scheme, 
even with the TVB limiter, shows a reduced accuracy for quite a large region around the 
shock. This indicates that, for a scheme which is globally oscillatory (like the central compact 
scheme), limiters can render it stable but may kill accuracy in smooth regions since there 
are oscillations there to suppress. On the other hand, the upwind compact scheme work 
well, with bigger errors for the TVD limiter near the smooth extremum which is close to the 
shock. The errors for the two dimensional case are similar and are not displayed. In the last 
plot, Figure 7, we show the surface of the two dimensional solution at t = 1 with 40 x 40 
points using the third order upwind method with TVB limiting. 


5 Concluding Remarks 

We have discussed a general framework to apply local limiters on compact schemes via the 
definition of a local mean. The resulting schemes are proven TVB (total variation bounded) 
in one dimension and maximum norm stable for multiple space dimensions. Numerical 
examples show that the base compact scheme should be upwind-biased in order to obtain 
high order accuracy after limiting for shocked problems. 
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Figure 1: L\ error versus number of grid points in log-log scale for smooth solutions. Stars: 
compact schemes without limiter; squares: with TVD limiter; diamonds: with TVB limiter. 
1(a): Third order upwind scheme, ID 
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1(b): Fourth order central scheme, ID 
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1(c): Third order upwind scheme, 2D 
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1(d): Fourth order central scheme, 2D 
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Figure 7: Surface of third order upwind compact scheme with TVB limiter for shocks, with 
40 x 40 points. 
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